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INTRODUCTION 

Solutions  to  problems  involving  systems  containing  nonlinear 
elements  are  seldom  clear-cut.   All  approaches  to  solutions  of 
problems  of  this  nature  are  laced  with  approximations  and  as- 
sumptions which,  while  simplifying  the  solutions,  are,  to  the 
degree  of  the  approximations  and  assumptions,  inaccurate. 

The  observation  (l)  may  be  made  that  the  coefficients  in 
equations  describing  physical  systems  are  never  known  to  a  high 
degree  of  accuracy.   These  coefficients  must  be  found  empiri- 
cally as  a  result  of  experimental  measurements  which  are  always 
subject  to  error.   Especially  in  nonlinear  systems  where  the  co- 
efficients are  functions  of  the  operating  conditions,  it  is 
difficult  to  determine  values  of  the  coefficients  to  a  high  de- 
gree of  accuracy.   Additionally,  the  physical  parameters  are 
often  subject  to  change  with  ambient  conditions  such  as  time, 
temperature,  and  the  like.   Changes  of  this  sort  are  probably 
not  included  in  equations  describing  the  system.   As  a  result, 
the  coefficients  are  subject  to  considerable  uncertainty.   Thus, 
depending  on  the  nature  of  the  nonlinearity  and  the  method  in- 
volved in  the  solution,  the  solution  may  or  may  not  be  of  a 
high  degree  of  accuracy.   Consequently,  there  is  reason  to  ques- 
tion whether  or  not  the  solution  finally  obtained  actually  ap- 
plies to  the  physical  system  under  study. 

The  obvious  question  which  arises  when  a  physical  system 
subject  to  the  above  observations  is  involved  is:   despite  the 
uncertainty  inherent  in  the  given  data  and  the  method  of  solution 


of  the  problems  involved,  is  there  an  approach  that  will  give 
enough  insight  into  the  problem  that  corrective  action  can  be 
taken  on  the  system  to  insure  a  reasonable  interpretation,  on 
the  part  of  the  system,  of  the  desired  performance? 

The  answer  to  this  question  lies  in  definition  of  the  word 
"reasonable".  '  In  general  the  minimum  requirement  for  a  reason- 
able performance  on  the  part  of  the  system  is  the  requirement 
that  the  system  be  stable.   This  in  turn  requires  definitions 
of  system  stability. 

Defining  Stability 

The  following  is  an  attempt  to  integrate  the  ideas  of 
Hughes  (2)  and  Cunningham  (l)  regarding  definitions  of  stability 
The  definitions  are  for  three  types  of  stability:   Asymptotic, 
orbital,  and  structural. 

Perhaps  the  original  essence  of  the  stability  idea  asks: 
if  a  system  is  initially  at  rest  or  else  operating  in  a  steady 
state,  either  of  which  is  to  say  that  the  system  is  operating  at 
some  equilibrium  point,  and  a  small  perturbation  is  applied  to 
the  system,  does  the  system  return  to  its  initial  state,  depart 
in  a  monotonic  increase  with  time  from  the  initial  state,  or 
achieve  some  ultimate  state  different  from  the  initial  state? 
At  first  glance  this  seems  to  be  a  fair  evaluation  of  what  can 
happen  to  a  system  and  for  a  linear  system  it  is  quite  adequate, 
the  key  words  being  "equilibrium  point". 


To  see  the  significance  of  this,  note  that  a  physical  sys- 
tem may  be  described  by  a  set  of  simultaneous  differential 
equations  of  the  form: 

dx-j/dt  =   x1   =   r1(x1,    x2,  x^,  .  .  .  ,  xn) 

dx2/dt  =  x2  =  f2(x-j_,  x2>    x3>  '     '     •>    -^n^ 


dx^/dt  =  xn  =  fn(x]_,  x2,  x^,  .  .  .,3^) 
with  t  the  independent  variable.   In  this  case  x^,  x2,  .  ..,  xn 
are  the  dependent  variables,  and  functions  f-^,    f?,  ...,  fn  are 
generally  nonlinear  functions  of  the  dependent  variables.   The 
simplest  equilibrium  points  are  those  points  where  x-^,  x2, 
. . . ,  x^   are  all  zero  simultaneously.   The  system  is  then  at 
rest  since  all  of  the  dependent  variables  are  constant  and  un- 
varying  with  time. 

For  a  linear  system,  functions  f-,,  f2,  ...,  fn  are  linear 
functions  only  and  when  the  derivatives  are  set  equal  to  zero; 
the  linear  functions  give  the  conditions  for  equilibrium  as: 
0  =  a11x1  +  a12x2  +  ai3x3  +  .  •  .,  ainxn 
0  =  a2]_x-^  +  a22x2  +  a2oXo  +  .  .  .,  a2nxri 


0  =   anlxl  +  an2x2  +  an3x3  +  •  '  "    annxn 
and  the  a.-  •  coefficients  are  constants  arising  from  the  param- 
eters of  the  physical  system.   Generally  the  determinant  of  the 
coefficients  does  not  vanish;  i.e., 
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In  this  linear  case,  then,  it  is  apparent  that  the  vari- 
ables Xj_  satisfying  the  system  equations  must  all  be  zero.   Thus 
a  linear  system  has  but  a  single"  equilibrium  point  at  which  all 
the  dependent  variables  vanish. 

For  a  linear  system  the  previous  definition  for  stability 
is  rigorous  since  the  necessity  for  considering  equilibrium 
points  other  than  zero  is  unnecessary.   For  a  linear  system  this 
type  of  stability  is  the  so-called  "asymptotic  stability". 

For  the  nonlinear  system,  f-,,  fp,  fo,  ••-,    fn  are  nonlinear 
functions  and  lead  to  nonlinear  equations  in  place  of  the  linear 
set.   These  nonlinear  equations  may  lead  to  solutions  for  the 
Xj_  which  are  not  zero  and  additionally,  more  than  a  single  set 
of  solutions  may  exist,  which  is  equivalent . to  saying  that  non- 
linear equations  may  have  many  equilibrium  points. 

From  the  above,  asymptotic  stability  which  characterizes 
the  system  as  stable  if  it  returns  to  the  initial  state,  un- 
stable if  it  continues  to  diverge  from  the  initial  state,  or 
"neutrally"  or  "conditionally"  stable  if  it  attains  some  new 
equilibrium  state,  does  not  describe  adequately  some  conditions 
which  may  appear  in  nonlinear  systems.   For  the  nonlinear  sys- 
tem it  is  necessary  to  specify  that  the  initial  disturbances  of 
the  Xj_  be  small  enough  to  keep  them  in  the  region  of  the 


equilibrium  point  in  question.   This  definition  may  be  shown 
graphically.   If  a  system  under  study  is  operating  about  some 
equilibrium  point  and  a  small  but  finite  perturbation  is  applied 
such  as  shown  on  the  phase  portrait  in  Fig.  1,  then  if  P  repre- 
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Fig.  1. 
sents  an  unperturbed  moving  point  on  the  phase  portrait  and  P' 
represents  the  corresponding  point  on  the  perturbed  portrait, 
and  if 
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the  system  is  said  to  be  asymptotically  stable.   The  motion 
shown  is  periodic  but  periodic  motion  is  not  necessary  for  the 
argument . 

The  significance  of  structural  stability  may  be  shown  as 
follows  (L(.)  :   for  nonlinear  problems  normally  encountered,  it 
is  often  necessary  to  express  some  nonlinear  function  of  a  var- 
iable and  its  derivatives  with  a  finite  number  of  terms  of  a 
power  series.   Consider,  for  example,  a  second  order  system 
whose  phase  portrait  is  described  by 

dVb   A(Va,  Vb) 


dVQ    B(V  ,  Vh) 


(1) 


with  A  and  B  simply  functions  of  Va  and  Vb.  Also  consider  a 
slightly  different  system  in  which  the  phase  portrait  is  de- 
scribed by 

dVb    A(Va,  Vb)  +  a(Va,  Vb) 


dV. 


B(Vfl,  Vb)  +  b(Va,  Vb) 


(2) 


with  a(VQ,  Vb)  and  b(Va,  Vb)  functions  of  Va  and  V-^  and  very 
much  smaller  than  A(Va,  Vb)  and  B(Va,  Vb)  ,  respectively.   If  the 
phase  portraits  of  these  two  equations  are  plotted  simultaneously 
with  the  same  initial  conditions,  the  plot  might  appear  as  shown 
in  Fig.  2. 
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Fig.  2. 

If  these  two  phase  portraits  are  qualitatively  similar  in 
nature,  the  original  system  may  be  defined  as  structurally 
stable.   The  above  differential  equations  are  really  expressions 
of  the  nonlinear  function  as  a  finite  polynomial  in  the  regions 
of  interest.   The  addition  of  the  terms  a(V  ,  V"b)  and  b(Va,  Vb) 
means  simply  that  the  coefficients  of  the  polynomials  used  to 
represent  A(V  ,  V\  )  and  B(V  ,  V^)  may  vary  slightly.   If  these 
slight  variations  markedly  change  the  system  behavior  as  evi- 
denced by  the  phase  portrait,  then  the  method  used  to  approxi- 
mate the  nonlinear  functions  is  of  doubtful  validity.   Use  of 


power  series  approximations  and  incremental  solutions  is  thus 
complicated  by  structural  instability.   The  indication  is  that, 
at  best,  incremental  solutions  must  be  much  more  precise  in 
order  to  get  a  reasonable ' determination  of  system  character- 
istics . 

Finally,  consider  a  system  having  a  steady-state  oscilla- 
tory motion  represented  as  a  closed  curve  in  the  phase  plane  as 
in  Fig.  3.   If  the  given  motion  is  slightly  perturbed  and  the 
new  motion  remains  within  the  immediate  vicinity  of  the  old 
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Fig.  3. 
motion,  then  the  system  is  said  to  be  orbitally  stable.   More 
rigorously  (2)  in  the  case  of  periodic  motion,  a  system  is  orbi- 
tally stable  if  for  a  given  perturbation  o  there  is  a  number  6 
such  that  the  distance  between  two  corresponding  points  on  these 
two  phase  portraits  is  never  greater  than  5.   In  view  of  this 
definition  and  the  definition  for  asymptotic  stability,  a  con- 
servative nonlinear  oscillating  system  is  not  orbitally  unstable, 
although  it  is  asymptotically  unstable. 


The  purpose  of  this  report  is  to  review  some  of  the  methods 
of  determining  stability  aspects  of  nonlinear  systems. 

THE  DESCRIBING  FUNCTION  APPROACH 

The  so-called  "describing  function"  is  an  application  of 
the  principle  of  harmonic  balance  to  various  nonlinear  elements. 
The  background  for  the  principle  of  harmonic  balance  and  the 
describing  function  is  discussed  in  Appendix  A.   Describing 
functions  for  many  types  of  nonlinearities  have  been  derived, 
as  in  Kuo  (3)  and  Gibson  (5).   As  an  example  of  the  derivation 
of  a  describing  function,  which  will  be  used  in  a  stability 
study,  take  the  case  of  a  transmission  element  saturating 
abruptly  ( 1) „ 

Many  physical  transmission  elements  have  the  property  of 
saturation.   The  output  quantity  is  related  linearly  to  the 
input  quantity  so  long  as  the  input  magnitude  is  less  than  some 
critical  value.   If  saturation  is  assumed  to  take  place  abruptly, 
the  input-output  relation  can  be  described  simply  in  terms  of 
several  linear  algebraic  equations.   If  the  input  instantaneous 
value  is  x  and  if  the  output  instantaneous  value  is  y,  they  may 
be  related  as  in  Fig.  ij.. 

These  relations  are 

-xc  <  x  <  +  xc    -   y  =  kx 
x  >  +  xc        -   y  =  +kxc 
x  <  -  xc        -   y  =  -kxc 
where  k  is  a  positive  constant.   Electronic  amplifiers  often 


saturate  abruptly  in  this  fashion 


Fig.  1^. 
The  describing  function  H  may  be  found  for  this  type  of 
element  by  assuming  that  the  input  quantity  varies  in  a  simple 
harmonic  fashion  as  x  ~  X  cos  cot,  where  X  is  the  amplitude  and 
co  is  the  angular  frequency.   The  output  quantity  will  also  vary 
in  a  simple  harmonic  fashion  if  X  <  xc,  in  which  case 
y  =  kx  cos  cot.   If  X  ?  xc,  the  output  quantity  is  not  simple 
harmonic  but  must  be  expressed  as  a  Fourier  series.   The  com- 
ponent of  fundamental  frequency  is  the  one  of  interest  and  it 
is  found  in  the  usual  manner.   If  the  input  and  output  quanti- 
ties are  plotted  as  functions  of  time,  the  result  is  shown  in 
Fig.  £. 
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Fig.  5- 
The  relationships  for  the  first  half  cycle  are: 


0  <   cot  <  -< 

°<    fr    COt    £    71    -    «< 
7t     -     <X    £    COt    <    71 


X  >  +  X, 


-xc  <  x  <  +  xx 


X, 


y  =  +kxc 

y  =  +kx  cos  cot 
y  =  -kxc 


X, 


where  »<  Is  the  angle  for  which  cos  «<  =  —  and  X  >  xc .   The 

x 
second  half  cycle  is  symmetrical  with  the  first,  and  thus  it  is 

unnecessary  to  consider  it  in  the  Fourier  analysis.   Because  of 

the  symmetry  only  cosine  components  appear  in  the  output  wave. 

The  amplitude  of  the  fundamental  cosine  component  is: 


Ai  = 


•7t-»< 


kx~  cos  cot  d(cot)  + 


'0 


kx  cos   cot  d(cot) 


'*< 


,7t 


(-kxc)  cos  cot  d(cot) 


(tc--0 


Carrying  out  the  integration, 
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cos  cot  d(cot)  -      cos  cot  d(cot) 

/(n-0 

+  kx        cos   cotd(cot) 
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2 
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-    sin  cot 

71 
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"l                                           1 
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+  -  kx 

—   sin  cot    cos  cot   +     •  cot 

71 
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</ 
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r    2 

kx  I  1  -  —  (»<  -  sin  »<  cos  «<) 

L      ^ 


If  X  =  x  ,  then  ^  =  0  and  A,  =  kx  .   If  X  >>  x  .  then 
c'  1     c  c' 


71 


-4  '  '  and  A-.  — >  

2       1      71 
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The  describing  function  is  then  the  ratio  H  =  —  and  is 

x 


X^  xc 
X  >  xn 


H  =  k 
H  =  k 


1  -  —  («*  -  sin  <a   cos  «<) 

71 


X, 


where  cos  =<  =  —  .   The  describing  function  depends  upon  the  in- 

x 
put  amplitude  but  not  upon  the  frequency.   It  is  a  real  number 

having  magnitude  but  zero  angle.   It  varies  as  a  function  of  the 

X 
ratio  —  as  shown  in  Fig.  6. 

xc 

The  describing  function  may  be  put  to  further  use  in  a 

stability  study  after  some  preliminary  remarks  about  Nyquist's 
criterion  and  the  stability  of  feedback  systems. 

Assuming  a  knowledge  of  the  Routh-Hurwitz  criterion  (7)  for 
stability  and  the  characteristic  equation,  this  criterion  may  be 
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Pig.  6. 
stated  simply  by  saying  that  the  requirement  for  stability  is 
that  the  characteristic  equation  and  its  coefficients  be  posi- 
tive, and  further  that  the  characteristic  equation  have  no  roots 
with  positive  real  par  us  (or  no  pure  imaginary  roots,  in  a  prac- 
tical sense).   The  Routh-Hurwitz  criterion,  however,  does  not 
give  any  information  concerning  methods  of  improving  the  system. 

The  Nyquist  criterion  possesses  the  following  features  which 
make  it  particularly  desirable  for  the  stability  analysis  of 
feedback  control  systems  (3). 

1.  It  provides  the  same  amount  of  information  on  the 
absolute  stability  of  a  feedback  system  as  the  Routh 
criterion. 

2.  In  addition  to  the  absolute  system  stability,  the  Ny- 
quist criterion  also  indicates  the  degree  of  stability 
of  a  stable  system  and  gives  information  about  how  the 
system  stability  may  be  improved  if  necessary. 

3.  The  Nyquist  locus  gives  information  concerning  the 
frequency  response  of  the  system. 

The  Nyquist  method  is  based  on  conformal  mapping  of  complex 
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quantities  and  stripped  of  its  rigorous  background  is  quite 
simple  to  use. 

Note  (1)  that  the  characteristic  equation  may  be  written 


n-1 


n-2 


an-l^  +  an  =  ° 


fU)  =  a0\n  +   a1\L1~±   +   a2AJ 
The  value  of  X  in  the  characteristic  equation  is  generally  a 
complex  number  and  can  be  represented  on  the  complex  plane  as 
in  Fig.  7  with  ~k   =   5  +  jco.   Any  value  of  X    that  is  a  root  of  the 
characteristic  equation  and  leads  to  an  unstable  solution  has  a 


Fig.  7. 

positive  real  part  and  would  be  located  on  the  right  half  of  the 
complex  \   plane.   This  region  is  shown  shaded  in  Fig.  7.   Its 
boundaries  can  be  traced  out  by  starting  at  the  lower  end  of  the 
jco  axis  where  X   =   -jco  - — >   -jso  ,  moving  up  this  axis  until 
A  =  jco — >+  j  c-3,  turning  clockwise  through  a  right  angle,  and 
returning  to  the  starting  point  along  a  semicircle  of  very  large 
radius.   At  the  starting  point,  a  second  clockwise  right  angle 
turn  is  needed  to  begin  retracing  the  original  path.   The  shaded 
region  is  always  to  the  right  of  this  boundary  as  it  is  traced 
in  the  direction  indicated. 

The  algebraic  function  f (A)  is  involved  in  the  characteristic 


Ik 


equation.   Plot  this  function  on  the  f  plane  as  shown  below  as 
X    traces  out  the  boundary  of  the  figure  just  described  for 
f(X)    =   X     +   gX   +   h,  g  and  h  constants.  X    is  allowed  to  be  pure 
imaginary  and  takes  on  successive  values  between  -joo  and  +jo£>  . 
Corresponding  values  of  f(jco)  are  calculated  and  plotted.   This- 
leads  to  the  heavy  curves  shown  on  the  X   plane  (Fig.  7)  and  f 
plane  (Pig.  8) .   The  dotted  curve  of  the  f  plane  is  found  by 
tracing  the  solid  curve  to  the  point  co — v  +  c*=?  f    turning  clock- 
wise through  a  right  angle  and  returning  to  the  point  co — >   -0*3 
along  a  circular  path.   As  the  boundaries  in  either  the  X   plane 
or  f  plane  are  traced  in  the  directions  indicated,  the  shaded 
areas  correspond  and  are  located  to  the  right  of  the  path. 


co 
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f  plane 


f(H=o 
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Pig.  8. 
Since  the  polynomial  f{X)    is  of  degree  n  in  X,    each  point 
on  the  f  plane  has  corresponding  to  it  n  points,  generally  all 
different,  in  the  X   plane.   By  the  same  token,  a  given  point  on 
the  plane  has  only  one  point  corresponding  to  it  on  the  f  plane 
Any  point  in  the  shaded  region  of  the  X   plane  must  have  corre- 
sponding to  it  a  point  in  the  shaded  region  of  the  f  plane. 
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These  points  represent  values  of  A.  which  would  lead  to  unstable 
solutions . 

The  characteristic  equation  is  f(X)  =   0,  which  is  repre- 
sented by  the  origin  of  the  f  plane.   Thus  the  roots  of  the 
characteristic  equation,  which  determines  properties  of  the 
solution  for  the  original  differential  equation,  are  represented 
by  those  points  in  the  X  plane  which  correspond  to  the  origin  of 
the  f  plane.   Instability  is  indicated  and  at  least  one  of  the 
roots  has  a  positive  real  part  if  the  origin  of  the  f  plane  is 
in  the  shaded  area. 

In  the  analysis  of  feedback  systems  the  equation  being 

PiU) 

studied  often  involves  fractions  of  the  form  f(X)  =  1  +  

P2U) 

where  P]_(X)  and  P2(X)  are  polynomials  in  X.   Here  f(X)  may 

become  infinite  for  certain  values  of  X  =  jco  such  that  Pg^) 

=  P2(jco)  =  0.   These  are  the  poles  of  f(X).   The  procedure  for 

mapping  is  to  avoid  the  poles  of  f(X)  by  following  a  small  semi- 

k 

circle  around  them.   If  f(X)  =  1  +  ,  where  k  and  2  are 

X(\   +  I) 
constants,  a  pole  is  located  at  X  =  0.   This  is  analogous  to 

the  case  already  discussed.   The  path  to  be  followed  is  shown 
in  Fig.  9a,  where  a  small  detour  is  needed  around  the  origin 
of  the  X  plane.   The  corresponding  path  in  the  f  plane  is  shown 
in  Pig.  9b,  where  a  large  semicircle  appears  corresponding  to 
the  small  semicircle  about  the  origin  of  the  f  plane.   In  each 
case  as  co  increases  from  negative  to  positive  values,  a  clock- 
wise right-angled  turn  is  made  as  the  semicircular  path  is 
entered.   Figure  9  again  indicates  a  stable  system.   Typically 
the  Nyquist  method  involves  considerable  numerical  calculation. 
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Another  remark  is. in  order  concerning  the  application  of 
the  Nyquist  criterion.   Block  diagrams  as  shown  in  Fig.  10  are 
frequently  used  in  describing  feedback  and  other  systems. 


*b 


*r\ — U~ 


H 


Fig.  10. 
The  system  above  may  be  represented  by  the  equations 

x2 


xl  ~  *0  "  x3  ~ 


0(D) 


=  Xq  -  H(D)x2 


G(D)xq 


or 


x2 


1  +  G(D)H(D) 


where  G  and  H  are  generally  functions  of  the  derivative 

d 
operator  D  =  —  . 

dt 

Asymptotic  stability  of  the  system  is  determined  by  what 
happens  in  the  absence  of  any  input  signal  to  the  system,  so 
that  Xq  =  0.   With  Xq  =  0,  x2(l  +  G(D)H(D))  =  0.   Typically,  a 
solution  for  this  differential  equation  is  assumed  as  x2  =  Xe   , 
where  X  is  an  arbitrary  constant  and  A  is  the  characteristic 
exponent.   Substituting  the  solution  into  (l  +  G(D)H(D))x2  =  0 
yields  f(X)    =   1   +   G(A)H(A.)  =  0.   Thus  the  system  stability  may 
be  tested  by  the  Nyquist  method.   In  applying  the  Nyquist  plot 
it  is  slightly  simpler  to  plot  not  f(jco)  but  rather  the  curve 
representing  just  the  product  G(  jco)H(  jco)  .   Stability  is  governed 
by  the  relation  of  this  curve  to  the  point  at  -1  +  jO .   In 
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addition,  any  physical  system  is  a  low-pass  device  so  that  the 
product  G(  jco)H(  jco)  approaches  zero  as  go  approaches  infinity. 
Therefore  the  curve  representing  the  product  closes  on  itself. 
The  system  is  stable  if  the  curve  for  G(jco)H(jco)  does  not  en- 
close the  point  -1  +  jO .   The  system  is  unstable  if  the  point 
is  enclosed. 

In  applying  the  Nyquist  criterion  as  described,  functions 
G  and  H  for  the  parts  of  the  system  become  merely  the  transfer 
functions  defined  with  simple  harmonic  variations.   If  some 
part  of  the  system  is  slightly  nonlinear,  its  transfer  function 
becomes  the  describing  function  as  previously  noted.   Provided 
that  the  nonlinearity  is  not  too  great  and  that  the  wave  forms 
in  the  system  are  essentially  sinusoidal  in  shape,  a  prediction 
of  this  kind  may  be  essentially  correct.   Testing  the  stability 
in  this  manner  is  open  to  question  if  the  system  is  such  that 
the  wave  forms  depart  considerably  from  a  sinusoidal  shape. 

In  view  of  the  preceding  discussions,  an  analysis  may  be 
made  of  a  system  containing  a  nonlinear  element  utilizing  the 
techniques  described  (1).   The  circuit  shown  in  the  block  dia- 
gram in  Fig.  11  consists  of  an  electronic  amplifier. with  a 
phase-shifting  network  connected  between  its  output  and  input 
terminals.   Since  any  practical  amplifier  saturates  if  the  mag- 
nitude of  the  signal  voltage  applied  to  it  becomes  too  large, 
the  amplifier  becomes  nonlinear  in  the  large  signal  node. 

Assume  at  first  that  the  amplifier  is  linear  without 
saturation  effects,  that  is,  operating  in  the  small  signal  mode. 
The  limiter  is  thus  not  considered,  and  therefore  e2  and  e^  are 
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Fig.  11. 
identical.   A  number  of  phase-shifting  networks  may  be  employed 
with  oscillators  of  this  type.   The  one  shown  is  sometimes  used. 
The  transfer  function  for  this  phase-shifting  network  is  (l\.)  : 

Gk                                          1 
H2(D)  =  3  =  = = 

e3    (RCD)-^  +  5(RCD)^  +  6(RCD)  +  1 

Attenuation  and  phase  shift  both  occur  on  a  signal  passing 

through  this  network  and  both  increase  as  frequency  increases. 

Specifically,  take  the  amplifier  as  a  single-stage  vacuum- 

e2 
tube* circuit  with  voltage  amplification  G  -  —  =  3?.   Typically, 

el 
this  amplifier  would  reverse  the  polarity  of  the  signal, 

justifying  the  negative  sign  on  eu    and  the  algebraic  signs  of 

the  equations.   In  this  case,  then,  the  stability  of  the  system 

is  given  by 

1  +  GH2U)  =  0 

and  the  test  for  stability  consists  of,  in  part,  plotting 

35 


GH2(jco)  = 


[_-5(RCa))2  +  l]  +  j[-(RCoo)3  +  6(RCoo)] 
as  co  varies  from  -<*a   to  +  o£>    .   This  curve  is  shown  as  Fig.  12. 


i  Imaginary 
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Pig.  12. 
As  co  increases  in  the' positive  direction,  the  point  -1  +  jO  is 
always  to  the  right  of  the  curve  and  is  completely  circled  by 
the  curve.   The  system  is  asymptotically  unstable,  and  oscilla- 
tion will  occur  with  increasing  amplitude. 

A  practical  amplifier  always  saturates  if  the  applied  sig- 
nal becomes  too  large.   The  saturation  effect  can  be  considered 
by  including  a  limiter  following  the  ideal  amplifier  as  was 
shown  in  Fig.  11.   A  simple  limiter  was  the  object  of  one  of  the 
preceding  discussions  and  the  results  will  be  used  here.   The 
limiter  is  assumed  to  saturate  abruptly  when  the  magnitude  of 
the  input  voltage  e2  becomes  too  large;  that  is,  exceeds  a 
critical  value  eQ,  so  that  the  following  relationships  apply. 
-ec  <  e2  <  +  ec:  e3  =  e2 

e2  ^  ec:      .  e3  =  +ec 

e2  £  -ec:  e3  =  -ec 

These  are  the  relationships  of  the  describing  function  example 
with  k  replaced  by  unity.   With  a  sinusoidal  input  voltage  of 
amplitude  E2  applied  to  the  limiter,  the  describing  function 
for  it  has  been  shown  to  be 
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E2  <    ec: 

H   =    1 

E2  >    ec: 

H   = 

where 

*c 

=  —      and    k   = 
Eo 

1    . 

1  -  —  (-<  -  sin  -<  cos  «<) 

71 


The  stability  of  the  system  including  the  limiter  is 
accordingly  governed  by  the  relation 

1  +  GH(E2)H2(A)  =  0  (1) 

This  is  the  equation 

f(X)  =  1  +  H1(A)H2(A)   with  E±   =   GH(E2) 
where  H(E2)  is  the  describing  function  for  the  limiter.   The 
equation  can  be  put  into  more  convenient  form  for  study  by 
writing 

1 


=  -H(E2) 


(2) 


GH2U) 

which  is  the  equation  for  steady-state  oscillation.   For  this 
example  G  =  35,  H2(M  was  previously  given  with  D  replaced  by 
X,    and  H(E2)  is  given  by  the  final  figure  (Fig.  3)  with  k  =  1 
and  X  =  E2.   In  Fig.  13  are  plotted  two  curves  representing  the 
two  sides  of  equation  (2)  with  X   replaced  by  jco.   The  left  side 
of  the  equation  is  a  curve  with  certain  values  of  the  quantity 
RCgo  indicated  along  it.   The  right  side  of  the  equation  is 

simply  that  portion  of  the  negative  real  axis  lying  between  the 

2 

point  -1  +  jO  and  the  origin.   Certain  values  of  the  ratio  — 

ec 
are  indicated  along  this  line.   Equation  (2)  is  satisfied  when 

the  two  curves  intersect,  and  conditions  determined  by  the 

intersection  point  correspond  to  steady-state  conditions  in  the 

nonlinear  circuit.   Essentially  the  effect  is  this:   because  of 
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the  asymptotic  instability,  the  amplitude  of  oscillation  in- 
creases from  an  initial  small  value;  however,  the  limiter  re- 
duces the  effective  amplification  until  at  steady  state  the 
effective  amplification  plus  limiter  is  just  sufficient  to  over- 
come attenuation  in  the  phase-shift  network. 


n  Imaginary  axis 


RCco  =   2.8 


RCa)=   2.1+5 
H(E2)    =    .83 
E2/e      =    1.1+3 


+  .2 


>  »  Real  axis 


-  -.2 


Fig.  13. 
The  intersection  of  the  curves  above  with  the  circuit 
parameters  used  in  the  example  indicates  that  steady-state 

oscillation  should  occur  with  a  frequency  of  RCco  =  2.1+5  and  an 

E2 
amplitude  such  that  —  =  1.1+3.   The  flat-topped  wave  of  the 

■Sc- 
out put  as  shown  in  the  section  on  the  describing  function  should 

be  somewhat  similar  to  the  limiter  output  wave  form.   The 

limiter  output  is  a  sinusoid  with  a  peak  1.1+3  times  the  limiting 

level.   Since  the  phase-shifting  network  is  a  low-pass  filter, 

the  wave  form  at  the  input  of  the  amplifier  should  be  nearly 

sinusoidal.   If  this  is  the  case,  the  prediction  based  on  the 
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describing  function  for  the  nonlinear  element  should  be  fairly 
accurate.   If  amplification  G  were  increased,  a  larger  quanti- 
tative error  could  be  expected.   As  can  be  seen  from  the  output 
in  the  section  on  describing  functions,  the  greater  the  gain, 
G  (X  in  that  discussion) ,  the  more  nearly  the  output  approaches 
a  square  wave  and  the  less  applicable  the  describing  function 
method  becomes. 

ANALYSIS  BY  MEANS  OP  SINGULAR  POINTS 

The  tunnel  diode  is  an  example  of  a  two-terminal  resistance 
in  which  the  instantaneous  current  ir  and  instantaneous  voltage 
ep  are  related  by  a  curve  similar  to  that  of  Pig.  II4.  (1,  9). 
Other  means  are  possible  to  obtain  this  curve.   Utilizing  the 
ideas  in  Appendix  B  concerning  singular  points  and  considering 
the  nonlinear  element  as  being  composed  of  several  linear  ele- 
ments, a  stability  analysis  can  often  be  made. 


,^ir=f(er) 


214- 


The  functional  relationship  is  i   =  f(ep)  and  may  possibly 
be  found  by  experiment  or  expressed  in  analytic  form.   In  the 
central  region  the  variational  resistance  is  negative,  implying 
a  source  of  power.   For  the  tunnel-  diode,  the  region  between 
points  b  and  c  on  the  curve  represents  a  decrease  in  the  tunnel- 
ing effect  with  subsequent  decrease  in  current  as  the  valence 
band  in  one  region  of  the  diode  is  raised  to  a  position  opposite 
a  forbidden  band  in  the  other  region  by  the  increasing  voltage 

(9).   The  variational  resistance  at  any  point  on  the  curve  is 

1 

defined  by  r  =  

dJL. 


de- 


Consider  the  circuit  in  Fig.  15  for  an  analysis  by  means 
of  singular  points  where  the  components  other  than  the  box  are 
linear  and  a  constant  voltage  E  is  applied  to  the  circuit. 


0 
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•R 


iC 


E  -±r- 


Fig.  15- 
Let  ipj  be  the  battery  current  and  i   the  current  through 
the  nonlinear  element.   The  current  through  the  capacitor  is 
thus  ( ipj  -  ir)  .   The  loop  equations  are 

1 
C 


-er  =  - 


f(er) 


ipfdt 


2$ 


de. 


dt 


-  f(eT.) 


(1) 


C  L 


and 


Ldi 


E  =  RiR  + 


R 


dt 


+   [ (iR  -  ir)dt 


diR     1 

=  -  (E  -  RiR  -  er) 

dt     L 


(2) 


di 


R 


The  singular  or  equilibrium  points  exist  where  =  0  and 


de. 


dt 


dt 

=  0  simultaneously.   Thus  at  the  singular  points: 


E  =  RiRS  +  erS 
iRS  -  f(er)  =  irc 
The  fact  that  at  a  singularity  iRg  =  irg  indicates  that  the 
capacitor  current  is  zero. 

A    conventional  technique  can  be  used  to  determine  the 
singularities.   From  a  point  on  the  er  axis  corresponding  to 

steady  voltage  E,  a  load  line  may  be  erected  with  slope  deter- 

1 
mined  by  resistance  R.   The  slope  of  this  line  is  -  —  .   The 

R 
intersection  of  f(er)  with  this  load  line  determines  the  singu- 
lar points.   The  singularities  are  thus  dependent  not  only  upon 
the  nonlinear  characteristic  but  upon  the  applied  voltage  and 
loading  resistor  as  well.   Cases  of  special  interest  have  an 
intersection  in  the  negative  resistance  region.   In  Fig.  1  two 
cases  are  shown;  the  three  intersection  and  single  intersection 
cases.   There  is  a  basic  difference  in  operation  in  the  two 
cases,  the  first  case  being  a  switching  operation  and  the 
second  case  leading  to  an  oscillation. 
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Following  the  routo  described  in  Appendix  B,  near  a 
singularity  the  substitutions  can  be  made: 

i,R  =  iRS  +  i    and    eR  =  erS  +  e 
where  i  and  e  are  small  changes.   Equations  (1)  and  (2)  thus 

become : 


de     1 
dt     C 


1 
i  -  (-)e 

r 


di     1 

—  -  (-)  C-Ri  -  e) 

dt     L 


'  (3) 


(w 


or 


di 
de 


1 
(-)(-e  -  Ri) 

L 


1  r  1 
-)\  -(- 


(5) 


(-)  -(-)e  +  1 
C  I      r 


From  (5)  coefficients  may  be  written: 

1        1         1  R 

a  =  -(-)C,  b  =  -  ,  C  =  ,  d  =  

r        C         L  L 


and  the  characteristic  roots  are: 


Xl>    A2  =  -\  -  [—      +  'J   ± 
2    \  rC     L/ 


11 

_i_ 

R 

Arc 

L 

+ 


k 


LC, 


f  R 

i     r 


,1/2- 


(6) 


An  analysis  can  now  be  made  for  various  conditions  in  terms  of 
the  characteristic  roots  X-i    and  ^-2* 

Now  consider  the  straight-line  approximation  to  the  true 
characteristic  shown  in  Fig.  l6 .   The  figure  may  be  divided  into 
three  regions.   Taking  as  a  first  case  the  situation  where  E  and 
R  are  large  enough  so  that  three  intersections  occur,  the  varia- 
tional resistances  are  -r^,  +r2,and  +r->  in  regions  1,  2,    and  3, 
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respectively.   The  slope  of  the  load  line  is  -l/R  in  all  regions 
Now  the  types  of  singularities  in  each  region  may  be  examined. 
In  region  1  the  slope  of  the  load  line,  -l/R,  is  less  in  magni- 
tude than  the  slope  of  the  negative  resistance  characteristic, 
-l/r-,.   Thus  R  >  Pn,   Substitution  into  equation  (6)  indicates 


if   /i       R\     F/i       Rr     ij. 

Xt,    \o  =  ""  \  ~  \  —  +  _ ]±'(  —  +  —        +  —    (positive    number) 
2   I     ArC        L/        [VrC        L/  LC 


l/2 


and  the  characteristic  roots  are  real  and  of  opposite  sign.   The 
singularity  in  this  region  is  a  saddle  and  unstable. 


Region  No. 2 


Region  No.  1 


Region  No.  3 


Pig.  16. 
The  variational  resistances  are  +T2    and  +ro  in  regions  2 
and  3,  respectively.   Substitution   into  equation  (6)  indicates 
that  the  singularities  in  these  regions  are  stable  and  are 
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either  nodes  or  foci  depending  upon  the  relations  among  the  cir- 
cuit element  values.   The  circuit  used  under  these  conditions 
can  be  the  basis  for  certain  types  of  triggering  circuits  char- 
acterized as  having  two  stable  states  separated  by  an  unstable 
state,  the  triggering  action  coming  from  suitable  influences 
applied  from  outside. 

Several  solution  curves  for  small  variations  in  current  and 
voltage  are  shown  sketched  in  Fig.  l6 .   Since  the  types  of  singu- 
larities are  known,  the  only  additional  knowledge  necessary  is 

di  di 

that  —  =  0  along  i  =  -e/r,  the  resistance  load  line,  and  —  =co 

de  de 

along  i  =  e/r,  the  characteristic  of  the  nonlinear  resistance 

element,  plus  the  slopes  of  m-^  and  rri2  near  each  singularity. 

As  noted,  in  regions  2  and  3  the  singularities  are  stable  and 

are  now  assumed  to  be  nodes.   The  curves  represent  the  way  small 

variations  i  in  the  current  in  and  e  in  the  voltage  ep  change 

with  time  following  some  small  initial  disturbance  of  the  system. 

Within  any  one  region  the  system  is  assumed  .linear  so  i  and  e 

are  not  necessarily  limited  to  small  values.   If  the  system 

starts  at  some  initial  point  it  does  not  necessarily  come  to 

rest  at  the  equilibrium  point  nearest  the  initial-  point , as  may 

be  seen  by  observing  point  P  in  Fig.  l6;  but  the  system  will 

come  to  rest  at  one  of  the  two  stable  points. 

If  it  is  specifically  assumed  that  the  .circuit  is  in  the 

stable  state  characterized  by  the  singularity  in  region  2  and 

voltage  E  is  then  raised  by  an  amount  AE  sufficient  to  cause 

the  load  line  to  shift,  a  possible  situation  is  shown  in  Fig.  17. 

AE  is  actually  a  large  voltage  for  the  condition  shown.   The 
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Fig.  17. 
singularity  in  region  3  has  moved  to  position  3'.   However, 
singularities  associated  with  regions  1  and  2  have  moved  outside 
these  regions.   If  the  straight-line  characteristics  of  the  non- 
linear resistance  are  extended,  however,  "virtual"  singularities 
can  be  located  as  shown  in  Pig.  17.   Virtual  singularity  1'  is 
associated  with  region  1  even  though  it  is  actually  located  in 
region  2.   By  the  same  token,  2'  is  associated  with  region  2 
even  though  it  is  located  in  region  1.   The  singularities  are 
of  the  same  type  as  before  the  change  in  E  and  solution  curves 
are  sketched  as  before. 

In  the  case  where  only  one  intersection  occurs,  R  <   r-,  . 
Figure  18  indicates  this  situation. 

The  singularity  is  located  in  the  negative  resistance 
region.   Extending  the  portions  of  the  characteristic  with 
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Region  1 


Region  3 


Fig.  18. 
positive  slope  two  virtual  singularities  may  again  be  found, 
this  time  in  region  1.   If  it  is  possible  for  the  components  to 
be  chosen  such  that  r^  <  L/RC,  the  magnitude  of  the  resistance 
in  the  region  of  negative  slope  is  less  than  the  impedance  of 
the  parallel  LCR  circuit  at  its  resonant  frequency.   Again  sub- 
stituting in  equation  (6)  indicates  that  X^_    an<^  ^2  mus"k  have 
positive  real  parts  in  region  1.   Since  R  <   r^,  the  singularity 
cannot  be  a  saddle  and  must  be  a  node  or  focus.   If  it  is  a 
focus, 


kLC/  V 


> 


l< 


1    R 


and  oscillations  are  implied.   For  the  virtual  singularities  the 
argument  is:   r  >  0,  \-^   and  X2   have  negative  real  parts  from  (6), 
and  the  operation  is  stable. 
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In  Fig.  18  it  is  assumed  that  an  unstable  focus  exists  with 
singularity  1  while  stable  foci  exist  for  both  2'  and  3'.   The 
solution  curves  relating  instantaneous  values  of  ip  and  er  are 
spiral  curves,  spiraling  outward  from  singularity  1.   When  the 
curves  cross  into  regions  2  or  3,  they  become  spiral  curves  go- 
ing inward  toward  virtual  singularities  2'  and  3' •   An  important 
result  of  these  actions  is  that  eventually  a  closed  solution  is 
reached,  as  indicated  by  the  elliptical  path  in  Fig.  18.   This 
closed  curve  is  the  "limit  cycle"  and  appears  only  in  systems 
having  nonlinear  negative  resistance.   Quantitatively,  the  ex- 
istence of  the  limit  cycle  may  be  justified.   If  initially  the 
amplitude  is  very  small,  the  solution  curve  is  entirely  in  region 
1  and  the  amplitude  must  grow.   Alternatively,  if  an  extremely 
large  amplitude  were  to  exist  initially,  the  solution  curve 
would  be  mostly  in  regions  2  and  3,    where  the  amplitude  decays. 
Such  a  small  portion  of  the  curve  would  be  in  region  1,  with 
increasing  amplitude,  that  the  net  effect  would  be  decay.   There- 
fore initial  large  amplitude  decays,  initial  small  amplitude 
grows,  and  there  must  be  some  intermediate  amplitude  which 
neither  grows  nor  decays. 

Experiments  with  a  circuit  similar  in  nature  to  that  of 
Fig.  15  containing  a  tunnel  diode,  voltage  source,  magnetic 
coil,  and  resistance  in  a  configuration  such  as  Fig.  19,  indi- 
cate that  at  least  the  inductance  or  capacitance  of  Fig.  15 
must  be  considered  for  the  switching  or  oscillation  operations. 

By  means  of  Fig.  17,  an  attempt  can  be  made  to  show  why 
these  actual  components  or  else  parasitic  reactances  must  be 
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diode 


Fig.  19. 
utilized.   Figure  17  indicates  that  at  least  one  reactive  ele- 
ment is  essential  for  an  explanation  of  what  takes  place  during 
the  interval  of  triggering.   During  this  interval  the  path  fol- 
lowed must  be  the  point  relating  instantaneous  values  of  ir, 
and  ep  must  be  the  characteristic  of  the  nonlinear  resistance 
element.   If  no  reactances  were  present,  the  point  relating  in- 
stantaneous values  of  i-^    and  e   would  have  to  move  along  the 
straight-line  path  of  the  resistance  load  line.   Without  re- 
actance, ip  must  equal  ir  and  points  corresponding  to  these 
two  currents  must  coincide.   As  previously  noted,  this  is  not 
the  case  except  at  singularities  and  triggering  between  singu- 
larities could  not  occur.   Thus  a  difference  between  paths  could 
not  occur  except  through  a  component  of  voltage  drop  across  the 
series  inductor  or  because  of  a  current  through  capacitor  C. 
Possibly  only  one  of  the  reactive  elements  is  actually 
necessary,  in  which  case  a  single  first-order  equation  is  suf- 
ficient to  describe  the  circuit.   The  equation  describing  the 
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circuit  if  L  =  0  is 


de       Re 

C  —  +  r  +  — 

dt       eR 


where  r  is  -Ft,  +^2,    or  +n,  depending  upon  the  region  of  opera 
tion.   There  is  a  single  characteristic  root: 


1  - 


(r  +  R) 


rR 


In  region  1,  r  <   0,  and  R  >    |  r   so  that  X  >  0  and  the  system  is 
not  stable.   In  regions  2  and  3,    r  >  0  so  that  X  <   0  and  the 
system  is  stable.   Assuming  the  inductance  present  but  C  =  0, 

yields 

L  [  —  J  +  (R  +  r)i  =  0 
Vdt/ 

. (R  +  r) 
In  this  case,  X   =   -   .   In  region  1,  r  <  0,  and  R>  j  r[, 

L 

so  that  again  X  <•   0  and  the  system  is  stable.   In  regions  2  and 
3,    r  >   0  so  that  X  ^   0  and  the  system  is  stable.   This  result, 
however,  is  a  contradiction  of  the  previous  analyses  which  were 
experimental  in  nature, and  the  conclusion  seems  to  be  that  a 
capacitance  in  parallel  with  the  voltage-controlled  resistance 
is  necessary  in  obtaining  the  operating  characteristics  of  the 
circuit . 

Finally,  it  seems  apparent  from  the  solution  curves  that 
whatever  reactances  are  present  in  the  circuit  play  a  part  in 
the  time  to  trigger  from  one  stable  state  to  the  other.  Gen- 
erally it  is  desired  to  keep  this  time  at  a  minimum  and  it  is 
possible  that  parasitic  or  stray  reactances  are  sufficient  to 
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serve  the  purpose  without  the  insertion  of  actual  reactances 
of  larger  magnitude. 


DIRECT  METHOD  OF  LYAPUNOV  AND  VARIABLE  GRADIENT 
METHOD  FOR  STABILITY  ANALYSIS 


Stability  information  about  a  system  may  at  times  be  ob- 
tained without  actually  solving  the  differential  equations  de- 
scribing the  system  (5,  6,  8,  10).   For  linear  systems  the 
Routh-Hurwitz  criteria  provide  an  approach  of  this  nature.   The 
core  idea  of  the  direct,  or  second,  method  of  Lyapunov  is  that 
it  is  sometimes  possible  to  form  functions  of  the  system 
equation(s)  and  time  which  possess  certain  properties  useful  in 
the  analysis  of  the  system.   Chief  items  of  interest  are  several 
theorems  leading  to  a  "Lyapunov"  function  or  "V"  function,  the 
generation  of  this  function,  and  the  analysis  of  the  function. 

The  equations  exhibited  in  the  introduction  form  what  is 
sometimes  known  as  a  "canonic  form"  and  are  again  presented 
here . 

dx1/dt  =  X]_(x]_,  X2,  ..  .  . ,  xn) 
dXp  /dt  —  Xp  { x-j  ,  Xp ,  .  .  .  ,  x.  ) 


dxn/dt  =  Xn(x1,  x2,  .  .  .  ,  xu) 
The  first  theorem  attributable  to  Lyapunov  regarding  these  equa- 
tions may  be  stated. 

For  a  set  of  first-order  differential  equations  of  the  form, 
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x.   :  a.  lXp  ^2 ,  •  •  •  ,  ^\J  **•'*•»    <-  >         •  •  >    n 
and  such  that, 

1.  xi  =  Xi(0,  0,  .  .  .,  0)  =  0     i  =  1,  2,  .  .  .,  n 

2.  The  functions  X^  are  continuous  with  respect  to  all 
variables  x-  in  the  entire  state  space. 

Then  there  exists  a  real-valued  scalar  function 
V(x-,,  X2>  .  .  .,  x^)  with  the  properties: 

a.  V(x3_,  X2,  .  .  .,  xR)  is  continuous  and  has  continuous 
first  partial  derivatives. 

b.  V(x-j_,  X2,  .  .  .,  x  )  >  0  except  when  x^  =  0  for 
i  =  1,  .  .  .,  n;  (i.e.,  V  is  positive  definite). 

c.  V(0,  0,  .  .  . ,  0)  =  0. 

2r 


d.  V(x1,  x2,  .  .  .,  xJ-900  for  (2-  x±*j  ~> c 

dv    n.  /pv  d*-i\  n 

e.  V  =  —  =  £_   •  1  <  -  e   21  *i 

dt    1=1  V^xi    dt  /        1=1 

(i.e.,  V  is  negative  definite). 
Then  the  system  is  asymptotically  stable  in  the  large.   Asymp- 
totic stability  in  the  large  assures  that  for  any  real  finite 
initial  conditions  on  the  system    the  output  of  the  system 
will  approach  the  equilibrium  state  of  x  -  0,  x  =  0,  as  t — »cx=> 
Some  further  definitions  are  necessary  here.   The  function 
V  is  called  positive  definite  or  negative  definite  in  a  given 
region  about  the  origin  if  at  all  points  in  this  region  it  has 
the  same  size  (positive  or  negative),  and,  except  at  the  origin, 
is  nowhere  zero.   The  function  V  is  called  semidef inite  if  it 
has  the  same  sign  throughout  the  region  except  at  certain  points 
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at  which  it  is  zero;  it  must  be  zero  at  the  origin  as  well.  The 
function  V  is  called  indefinite  if  in  the  given  region  about  the 
origin  it  takes  on  varying  signs. 

As  a  simple  example  of  the  use  of  the  theorem,  consider 
the  V  function  (without  considering  for  the  moment  how  it  was 
generated)  : 

1     2  2 

V  =  —  (Ipt-,  +  x]_x2  +  x2^ 

It  may  be  readily  verified  that  conditions  (a),  (c),  and  (d)  of 
the  theorem  are  satisfied.   The  question  arises  then  as  to 
whether  the  quadratic  form  above  is  positive  definite.   For  the 
general  quadratic  form 

F(x-,,  Xp)  =  ax-,  ^  +  2bx-,Xp  +  cxp 

and  a  necessary  and  sufficient  condition  for  it  to  be  positive 

p 

definite  is  that  a  >  0  and  ac  -  b  >  0.   For  the  form  above, 

a  =  1.33,  b  =  -33,    a^id  c  =  .33,  and  condition  (b)  of  the  theorem 
is  satisfied.   In  order  to  determine  whether  condition  (e)  is 
satisfied : 

»  • 

8    .     X]_X2    xpx!    2 

V  =  —  X]_xj_  +  +  +  —  xpxp 

3         3      3     3 

In  general,  there  will  be  a  relationship  between  x^_  and  xp . 
Let  Xo  =  (~3x?  "  3x-,  )  ,  Xp  -  x-j.   Thus 


lx2    2 


p    xl^    -* 
V  =  -(x-,   +  +  -  xp  ) 

3  3 

The  last  line  above  is  also  a  quadratic  form  and,  in  addition, 
is  negative  definite.  Since  all  'the  conditions  of  the  theorem 
are  satisfied,  it  is  asymptotically  stable  in  the  large  by 
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Lyapunov's  direct  method. 

An  analogy  may  be  drawn  between  the  potential  energy  of  a 
system  and  the  Lyapunov  function,  V.   The  positive  V  and  nega- 
tive  V  correspond  to  a  system  that  dissipates  energy.   The  re- 
sponse to  any  initial  condition  will  cause  the  system  to  dissi- 
pate energy  until  the  energy  is  zero.   Conditions  (1)  and  (2) 
of  the  theorem  are  satisfied  by  this  condition.   Herein  lies  a 
possible  key  to  the  generation  of  V  functions  but  there  are  a 
number  of  analytical  methods  of  generating  these  functions 
which  appear  more  promising. 

A  geometric  interpretation  may  be  placed  on  the  theorem 
and  this  interpretation  can  be  most  easily  shown  in  two-dimen- 
sional space.   Assume  a  V  function  of  two  variables  that  satis- 
fies conditions  (a),  (b) ,  (c),  and.(d)  of  the  theorem.   If 
V(x,  ,  Xp)  is  set  equal  to  a  constant  C,  then  the  resulting  equa- 
tion describes  a  closed  bounded  curve  in  the  x-,,  Xp  (state) 
plane.   If  V(x-,,  Xp)  is  a  quadratic  form,  then  the  curves  are 
ellipses.   If  C  is  allowed  to  have  different  values,  i.e., 
0  <  C-j_  <  O2  <  Co,  ...,  V(x-]_,  X2)  would  appear  as  in  Fig.  20. 

If  the  time  derivative  of  V(x-,,  Xp)  is  taken,  the  result  is 
dv    PV    dxx    2V    dx2 
dt   ^x-l   dt     2x2   dt 

Consider  a  phase  trajectory  crossing  the  curve 
V(x-,,  Xp)  =  G,  at  point  P.   The  partial  derivatives  may  be 
written: 


38 


^^-Phase  trajectory-- 
Cl   stable  system 


-Normal  to  V=Co  curve 
at  P  and  positive  in 
direction  shown 


Thus 

dv 
dt 


Fig.  20. 
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fdl 


fdV 


2' 


+ 


^   LV^X1 
5V 


1/2 


cos 


'.V  2 


9*2 


Jx2/    - 

2W2 


f 


+ 


>li/2  r 


cos  0 


+ 


«?x2/ 


dx-i  dx2 

cos  $  +  cos  9 


dt 


dt 


The  second  factor  on  the  right  represents  a  projection  of  the 
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tangential  velocity  along  the  trajectory  at  point  P  onto  the 

normal  of  the  V(x-,,  Xp)  -  C,  curve  at  point  P.   It  is  obvious 

that  if  V  is  negative,  then 

dx-[_         dxp 

cos  0  +  cos  9 

dt  dt 

must  be  negative  and  that  the  point  P  crosses  the  V(x,,  Xp) 
curves  from  outside  to  inside.   This  is  sufficient  to  assure 
that  as  t — 5»  o° ,    the  state  variables  x-^    and  X2  will  approach 
zero  and  the  discussion  may  be  considered  a  rough  proof  of  the 
theorem  for  the  two-dimensional  case. 

If  a  system  is  unstable,  Lyapunov's  direct  method  may  be 
used  to  establish  this  also.   A  theorem  to  do  this  is  identical 
to  the  theorem  previously  stated  with  the  exception  that  in 
condition  (e),  V  is  required  to  be  positive • definite  (that  is, 
of  the  same  sign  as  V) .   This  theorem  may  be  stated  in  a  some- 
what different  manner  also.   For  a  system  described  by  the 
canonic  set  previously  shown,  if  there  exists  a  real  valued 
function  V(x-,,  Xp,  .  .  .,  x^)  with  the  following  properties, 

(1)  V(x-,,  Xp,  .  .  .,  x  )  is  continuous  and  (2)  the  time  deriva- 

dVx 
tive  —  is  negative  definite,  then: 

dt 

1.  The  system  is  unstable  in  the  finite  region  for  which 

V  is  not  positive  semidef inite . 

2.  The  response  of  the  system  is  unbounded  as  t — >  oo 
if  V  is  not  globally  positive  semidef inite . 

A  graphical  interpretation  could  be  made  of  this  also  as  for  the 
stability  theorem  with  the  net  effect  of  showing  the  point  P  of 
the  phase  trajectory  crossing  the  C.  curves  from  the  inside  out. 
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■     A    definite  point  should  be  made  concerning  the  stability 
criteria  obtained  by  Lyapunov's  direct  method  which  is  that  they 
are  sufficient  to  establish  stability;  they  are  not  necessary 
criteria.   This  is  much  the  same  as  saying  that  if  a  function 
cannot  be  discovered  which  will  either  satisfy  the  conditions 
for  stability  or  instability  does  not  mean  that  the  function 
does  not  exist;  further,  a  Lyapunov  function  which  establishes 
a  stability  or  instability  for  a  system  is  not  unique. 

It  is  apparent  from  the  stability  and  instability  theorems 
that  the  major  problem  is  tne  actual  generation  of  V  functions 
in  order  to  test  for  stability  or  instability.  A   promising  ap- 
proach to  the  generation  of  V  functions  is  called  the  variable 
gradient  method.   The  method  is  summarized  below. 

If  a  system  is  stable  in  a  given  state  space,  then  a  V 
function  exists  for  the  system.   This  statement  follows  the 
stability  theorem  and  discussion  following  it.   Provided  only 
that  V  exists,  its  gradient  VV  exists.   The  last  statement  is 
discussed  in  reference  (6)  and  a  discussion  of  the  proof  given. 
Thus  given  the  gradient,  both  V  and  V  may  be  calculated. 
Actually,  the  approach  is  round  about  as  would  be  suspected  from 
the  use  of  VV  rather  than  V.   Without,  as  yet,  specific  know- 
ledge of  "V  V  except  that  it  will  be  some  function  of  the  state 
variables,  we  proceed  to  determine  V  from  \?V.      Thus 

<?V   dx1        dV        dx2         3V        dxn 

V  =  +  +    .     .     .    (1) 

^x-^dt  <?x2dt  <9xndt 

The    system   canonic    equations   may  now   be    inserted   for   the    x.j_    and 

V   put    into   vector   form: 


kl 


V  =  V  V  •  x  (2) 

Now  V  may  be  obtained  by  the  line  integration  of  equation  (2). 

r(x1,X2, • • . ,xu) 


V  = 


V  .  dx 


(3) 


'0 


This  line  integration  is  independent  of  the  path.   The  sim- 
plest path  of  integration  and  perhaps  the  easiest  for  actual 
integration  is: 

/xl  fx2 

V  =±   I   V1(yi,  0,0,  ...  0)dV1   +     I   V2(xx,  /2,  0,0,  ...0)d/, 


x 


n 


+  .  .  .  + 


Vn(x-L,  x2,  ...,  xn-1,  /n)d/t 


n 


(W 


o 


where  the  subscripts  on  the  V^  under  the  integral  sign  refer  to 
the  rows  in  the  matrix  representation  of  V  V.   Reference  (10) 
discusses  a  uniqueness  theorem  which  states  that  for  a  unique 
scalar  function  V  to  be  obtained  by  a  line  integration  of  a 
vector  function  VV,  the  generalized  curl  equations 


or 


3  x  j     3  *i 


x  vv  =  ° 


(5) 


must  hold.   This  is  the  n  dimensional  representation  of  Stokes' 
theorem.   Proceeding,  an  arbitrary  column  vector  VV  is  con- 
structed with  coefficients  which  are  themselves  functions  of 
the  state  variables. 


/'allxl  +  a12x2  +  •••  +  al 


nxn  ) 


VV  =   { 


21  1     22  2  2n  n 


(6) 


anlxl  +  an2x2  +  • • •  +  annx 


n 


k2 


Reference  (6)  discusses  the  coefficients  a.,  which  are  restricted 
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to  functions  of  x.  only  while  ann  is  set  equal  to  2  to  simplify 
the  process  of  closedness.   The  a •  .  are  composed  of  a  constant 
part  ajik  and  variable  part  a  i  ■?  v  ( *]_ ,  x2,  ...,  *n-±)  • 

In  summary,  the  outline  for  formal  application  of  the 
method  is: 

1.  Assume  a  gradient 'of  the  form  of  equation  (6). 
,  2.  Prom  the  variable  gradient  form:   V  =VV  •  x 
( equation  2)  . 
3.  In  conjunction  with  and  subject  to  the  requirements 
of  the  generalized  curl  equations,  equation  (5), 
constrain  dV/dt  to  be  at  least  negative  semidef inite . 
[|_.  From  the  known  gradient,  determine  V  and  the  region 

of  closedness  of  V. 
5.  Invoke  the  necessary  theorem  to  establish  stability. 
As  an  example,  following  the  steps  indicated  above, 
consider  the  system  shown  in  Fig.  21.   Let  the  nonlinear 


r(t)=p 


t  >  0 


^ 

NL 

3 

G 

-x 

lement  be  y  =  x^  and  G  = 


Fig.  21. 
1 


If  it  is  assumed  that 


S(S  +  1) 
the  equations  of  motion  written  in  state  variable  form  are 


w 


Xn  = 


1  =   x2 

Xp  ~  —  Xp  —  X-i 


then: 
Step  1 


VV  = 


allxl  +  fl12x2 


8  p  -i  X-i   t  8ppXp 


"W- 


Vv, 


Step  2 


V  =  (a-nXn  +  a-,pXp)x-i  +  (ap-,x-i  +  2x9)x 


11A1    a12^2y^l 


21^1 


^2'  2 


V  =  Xnx2^all  "  2xl   ~  a21^  +  x2  ^a12  ~  ^^  ~  a21xl 
Step  3 

To  constrain  V  to  be  at  least  negative  semidef inite,  the 

coefficient  of  x-.  x?  in  the  last  equation  in  Step  2  can  be  set 

equal  to  zero  and  also  0  <  a-^  ^  2.   a2]_  can  be  any  positive 

numb  e  r .   Thu  s 


V  =  -x22(2  -  a12)  -  ai2xi  >      an  =  a21 


+  2x- 


and 


Vv 


a21xl  +  2xl   +  a12x2 


S  q  -i  ^v-]    "t"   ^  -X-Q 


0  =  a12  <  2 


For  simplicity  let  a-,p  =  1 
Then 

2  x2       ^  xl 


implies 


kh 


1  = 


P(a21x1) 


a21x  +  a-21v  +  xl 


{a21v} 


and  if  a21v  ~  ^  an<^  a21x  ~  ^'  ^e  a^ove  relationship  is 
satisfied.   Finally, 


W 

Step  ij. 

V  : 


2x-^  +  x-]_  +  X2 

X.   +  2X0 


/Xl 

=  f   (2Xl3 


x2 


+  x-^)dx^  +j    (x-i  +  2x2)dx2 


0 


3 


xl^   xl 

V  =  +  +  xlx2  +  x2 

2      2  . 

Step  5 

Since  V  is  positive  definite  and  V  is  negative  definite, 
the  conclusion  is  that  the  system  is  globally  asymptotically 
stable.   Other  valid  choices  of  a-,  2  yield  other  valid  V  func- 
tions . 

In  summary,  autonomous  systems  seem  to  be  the  largest  field 
of  application  of  the  variable  gradient  technique.   The  method 
is  applicable  to  single  valued  continuous  nonlinearities  where 
the  nonlinearity  is  known  as  a  polynomial,,  a  specific  function 
of  x  or  a  curve  determined  from  experimental  results.   The 
method  generates  V  functions  to  suit  the  problem  at  hand. 
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APPENDIX  A 
The  Method  of  Krylov  and  Bogoliubov 

1.  General.   The  method  of  Krylov  and  Bogoliubov  (1,  2,    $) 
is  a  series  approximation  technique  for  determining  the  free 
periodic  oscillations  of  second  order  (and  with  a  more  general 
approach,  higher  order)  systems.   It  is  of  some  interest  not 
only  as  an  introduction  to  their  principle  of  harmonic  balance 
and  the  describing  function  technique  which  is.  based  upon  it, 
but  also  for  itself. 

Since  the  method  in  its  more  general  application  is  an 
equivalent  linearization  process  leading  in  some  cases  to  a 
kind  of  "transfer  function"  for  a  nonlinear  system,  it  will  be 
of  some  merit  to  briefly  review  the  linear  transfer  function 
concept. 

The  transfer  function  is  defined  when  simple  harmonic  var- 
iations with  respect  to  time  exist.   It  is  the  complex  ratio  of 
the  resulting  response  to  a  driving  force.   Depending  upon  the 
particular  situation,  the  driving  force  and  response   may  be 
either  currents  or  voltages  in  any  combination  (or  other  physical 
analogs  to  these  quantities) .   The  transfer  function  may  there- 
fore have  dimensions  of  impedance,  admittance  or  pure  numeric 
and  consists  of  real  and  imaginary  parts,  or  magnitude  and 
angle.   Both  parts  are  generally  functions  of  frequency  but  not 
of  amplitude  for  a  linear  system.   Transfer  functions  are  usually 
not  difficult  to  derive;  for  example: 
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where   E-^  =  input  voltage 

E2  =  output  voltage 

co   =  frequency  of  input 

m   =  jco 
In  the  course  of  analyzing  a  nonlinear  system,  it  is  often 
desirable  to  separate  the  linear  and  nonlinear  elements  so  far 
as  possible.   The  transfer  function  may  be  found  in  the  usual 
way  for  linear  elements.   A  kind  of  equivalent  transfer  function 
is  found  for  the  nonlinear  elements  considering  fundamental  com- 
ponents only.   This  equivalent  transfer  function  is  called  the 
describing  function  for  the  element.   Generally  it  is  a  function 
dependent  upon  amplitude  which  may  or  may  not  depend  upon 
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frequency .       It  leads  to  results  of  useful  accuracy  and  has  mean- 
ing only  in  the  event  that  the  nonlinearity  is  such  that  the  re- 
sponse to  a  simple  harmonic  driving  force  does  not  itself  dif- 
fer too  much  from  being  simple  harmonic. 

2.  The  Krylov  and  Bogoliubov  Method  for  Second-order  Systems 
Consider  the  following  equation: 

x  +  P(x,  x)  =  0  (1) 

x  =  the  control  (state  variable 

F(x,  x)  =  a  nonlinear  function  in  which  linear  terms 
are  contained  if  they  arise. 
Put  F(x,  x)  in  the  following  form: 

F(x,  x)  =  u)02x  +  u.f(x,  x)  (2) 

Equation  (l)  may  now  be  written: 

x  +  co02x  +  u-f(x,  x)  =  0  (3) 

The  approximations  below  require  that  u.  be  small.  A    solu- 
tion is  desired  of  the  form: 

x  =  a  sin(a)Q,t  +  <i) 
x  =  a  coq  c  o  s  (  coq  t  +  (f)) 
with 

a  =  a(t)  (6))   unknown 

0  =  0(t)       .  (7)/ 

For  equation  (5)  to  hold,  special  conditions  must  be  im- 
posed.  For  [i   =   0  in  (3),  the  solution  is  simple  harmonic  motion 
of  the  form  given  in  (1+)  with  a    and  <j)_   constant.   For  u.  ^  0,  the 
solution  is  no  longer  simple  harmonic  motion.   Functions  (6) 
and  (7)  permit  correction  of  nominal  frequency  o)q  and  if  they 
can  be  found,  equation  (I4.)  can  be  rigorously  correct.   This  is 
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no  mean  task  and  an  approximation  to  a(t)  and  $(t)  is  usually 
sought. 

Apparently  the  actual  frequency  at  any  time  is: 

0)  =  co0  +  —  w       (8) 

dt 

What  conditions  are  necessary  if  (5)  is  to  hold?   First 

differentiate  equation  (I4.)  with  respect  to  time. 

da  '  d 

x  =  —  sin(o)0t  +  0)  +  a  /cos(o)0t  +  6)  (coq  +  — 

dt  L  '        dt  ' 

da  d0 

x  =  —  sin(a)Qt  +  (j))    +   a  —  cos(o>Qt.  +  0)  +  ao)oCos(o)Qt  +  ffi)   (9) 

dt  dt 

If  (5)  is  to  hold: 

da  d<b 

—  sin(a)0t  +  $)+a  —  cos(o)0t  +  ®)=0  (10) 

dt  dt  ' 

Also  differentiate  (5)  with  respect  to  time. 

:.  da  d<P 

x  =  —  coq    cos(coQt   +  0)    -    a   Ua  j  sin(o)Qt   +   0)  (coq    +   ■ — ) 

dt  ^  dt 

da 


x  =  —  coq    cos(a)ot   +  0)    -    a   o)q      sin(o)Qt   +  0) 


-    a   goq  —   sin(o)Qt   +  0)  (11) 


dt 
Insert    (10)    and    (11)    into    (3)    with   z   =    (cont  +  0),    yielding: 


■'0 


da  .  o  d^ 


2 
dt   ."  "   dt 


o)q    cosz    -    acoQ      smz    -    ao)Q   —    sinz   +   o)q      a    smz 


+  u-f(a  sinz,  ao)Q  cosz)  =  0 

da  d(j) 

—  0)0  cosz  -  ao)Q  —  sinz  =  -u.f(a  sinz,  ao)Q  cosz)        (11a) 

dt  dt 
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da  d(f 

—  sinz  -  a  —  cosz  =  0  (10) 
dt          dt 

Multiply  (10)  by  o)q  cosz  and  (11a)  by  sinz  and  add, 
yielding: 

d(j)  1 

—  =  |xf(a    sinz,    scoq    cosz)sinz  (12) 

dt        aa)Q 

Multiply  (10)  by  goq  sinz  and  (lla)  by  cosz  and  add, 

yielding: 

da     1 


dt     (jOq 


|af(a  sinz,  ao)Q  cosz)cosz  (13) 


a(t)  and  $(t)  could  not  be  found  by  integration,  but  an  approxi- 


mate solution  for  the  set  must  usually  suffice.   Equations  (12) 

da  d(j) 
and  (13)  indicate  that  if  u-  is  small,  —  and  —  are  small,  and 

dt  dt 
that  a(t)  and  ^(t)  are  slowly  varying  functions  which  may  be 

considered  to  be  constant  at  their  average  values  over  a 

single  cycle.   Thus: 

da   _     1  f2% 

—  =  I  u-f(a    sinz,    a6)Q    cosz)cosz    dz  ( ll\.) 

dt        2-jigoq     /q 

dp  _        1  f2li 

—  == /  u-f(a    sinz,    ao)Q    cosz)  sinz    dz   .  (15>) 

dt         2naco0    /0 

Under   the    integral,    _a    is    assumed   constant    so    the    integra- 
tion   is   not    overly   difficult.      If    (2)    is    integrated: 
F(x,    x)     =    GOq    x   +    u-f(x,     x) 
F(x,    x)    -    C0q    x  =   u-f (x,    x) 
P(a    sinz,    aooQ    cosz)    -   coq    ( a    sinz)    =   u.f(a    sinz,    acoo    cosz) 
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da 

dt 


2ti 


27lO)0      ) 


P(a    sinz,    aa)Q    cosz)cosz    dz 
2         /2it 


0)0 


"O271    /0 


a    sinz    cosz    dz 


da 

dt 


2% 


=    -   I        P(a    sinz,    ao)Q    cosz)  cosz    dz 

27lO)0     /0 

Employing   the    square    of    (8)    with    (15)    yields: 


(16) 


9      .       o  d0         /d0^ 

0)^   =  o)0      +    2oo0   —  +      — 

dt        V  dt, 


2ti 


co^    =   goq      +    2o)q 


2Ttao)o 


P(a    sinz,    ao)Q    cosz)dz 


'0 


..    2     ,.271 


na       j 


a    sin2z    dz 


0 


a)2   = 


1 


2n 


F(a    sinz,    aa)Q    cosz)cosz    dz 


(17) 


Tta  /0 

With  (l6)  and  (17)  the  response  of  a  system  of  the  form  of 
(1)  may  be  found  in  the  form  of 

x  =  a ( t) sinz ( t) 
where  a(t)  and  z(t)  are  given  by 
da      a 
dt 


2co 


b(a,  co0) 


(18) 


0 


dz 


dt 


=  go2  =  g(a,  oo0) 


(19) 


and  where,  approximately, 


'"271 

F(a    sinz,    acoQ    cosz)cosz    dz 
Tia    lQ 


b(a,    o)n)    =  — 


(20) 
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l  f^ 
g(a,  coq)  =  —     F(a  sinz,  acog  cosz)sinz  dz  (21) 

na  /0 

3.  Generalization  of  the  Method.   In  the  above  form,  solu- 
tion may  be  obtained  for  nonlinear  equations  in  which  the  non- 
linearity  was  of  the  form 

F(x,  x)  =  g(a,  co)x  +  b(a,  co)x 
and  as  developed  applied  to  second-order  systems.   It  is  appar- 
ently possible  to'  obtain  solutions  of  a  similar  nature  to  higher 
order  systems  with  a  more  general  approach;  however,  for  demon- 
stration of  the  method  the  second-order  system  is  illustrative. 

Now  consider  the  nonlinear  function  y  =  F(x,  x)  and  let 

x  =  a  sin  cot .   Then  : 

y  =  F(a  sin  cot,  aco  cos  cot) 

Let  u  =  cot  and  expand  in  a  Fourier  series. 

1  f2^ 
y  =  —    F(a  sin  u,  aco  cos  u)du 

f1  (2%  ,  \ 

+  j—   J       F(a    sin  u,    aco    cos   u)sin  u   dur    sin  cot 

U  'o  ) 

/i  (2%  \ 

+  \  —  J      F(a  sin  u,  aco  cos  u)cos  u  dur  cos  cot 

+  higher  harmonics 
If  the  nonlihearity  is  symmetrical,  the  first  integral  in  the 
series  is  zero.   This  is  not  always  so.   For  a  symmetrical  non- 
linearity,  neglecting  higher  harmonics,  a  function  is  achieved 
of  the  form: 

b(  a,  co)  x 

y  =  g(a,  co)x  +  (22) 

co 


$k 


where  x  =  a  sin  cot,  x  =  aco  cos  cot  from  the  last  section;g(a,  go) 
mid  b(a,  co)  are  given  by  equations  (21)  and  (20),  respectively. 

Equation  (22)  is  linear  even  though  the  original  y  =  F(x,x) 
is  nonlinear.   In  connection  with  this,  note  that  what  has  been 
found  thus  far  by  the  method  of  Krylov  and  Bogoliubox  are  ap- 
proximation functions  which  replace  a  nonlinear  function  by  an 
equivalent  linear  function,  and  thus  the  appellation  "equiva- 
lent linearization".   The  approximations  are  useful  if  the  non- 
linearity  [i  is  small  so  that  the  amplitude  and  phase  charac- 
teristics of  the  nonlinear  element  may  be  considered  constant 
for  some  period  of  time.   In  other  words,  a  nonlinear  element 
may  be  approximated  by  an  equivalent  linear  one  and  the  approxi- 
mation is  useful  if  the  original  nonlinearity  is  not  too  large. 
The  approximation  may  be  accomplished  by  means  of  a  truncated 
Fourier  series. 

l\..    Example .   The  object  of  equivalent  linearization  in 
conjunction  with  the  principle  of  harmonic  balance,  then,  is  to 
choose  the  linear  replacement  element  such  that  the  fundamental 
slnz  and  cosz  components  are  the  same  for  both  linear  and  non- 
linear elements  under  a  simple  harmonic  motion.   The  equivalent 
linear  element  may  be  denoted  by  a  function  which  is  a  "describ- 
ing function"  for  the  nonlinear  element  subject  to  the  restric- 
tions of  the  approximations. 

As  an  example,  consider  a  combination  of  electrical  dis- 
sipatlve  and  reactive  nonlinear  elements  (or  possibly  just  a 
term  in  a  differential  equation)  for  which  the  voltage  is  some 
function  of  both  the  current  and  its  first  derivative.   Thus: 
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di 
v  -  f(i,  — ) 

dt 

For  an  equivalent  linearization  of  this  function,  the  linear 

function  is: 

di 
v  =  gi  +  b  — 
df 

If  it  is  assumed  that  the  current  through  the  element  is  sinu- 
soidal and  described  by 

i  =  I  c  o  s  ( w  t  +  0)  =  I  c  o  s  z 
the  voltage  across  the  nonlinear  element  or  term  is: 

v  =  fl  I  cosz,  -col  sinzj 

The  fundamental  sine  and  cosine  terms  are  respectively,  . 

1  (2* 
Vs;l  =  —  )   f(l  cos  z,   -col  sin  z)  sin  z  dz 
n  /0 

1  f2* 
Vc-|_  =  —  I   f(l  cos  z,  -col  sin  z)cos  z  dz 

The  voltage  across  the  equivalent  linear  element  is: 

v  =  (gi  cos  z)  -  (bco  I  sinz) 
By  the  principle  of  harmonic  balance  the  values  of  g  and  b 
may  be  found. 

.   1  z2* 

g  =  —  J   f(l  cos  z,  -co  I  sin  z)cos  z  dz 

1  /2* 
b  =  —  I   f(I  cos  z,  -co  I  sin  z)sin  z  dz 

di 
The  function  v  =  gi  +  b  —  is  thus  a  function  describing  the 

dt 
relationship  between  current  and  voltage  for  this  combination  of 

elements,  or  a  "sort  of"  transfer  function. 
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APPENDIX  B 
Analysis  of  Singular  Points 

The  analysis  of  the  singular  points  of  a  differential  equa- 
tion is  an  extension  of  phase  plane  analysis  and  can  be  useful 
in  determining  the  properties  of  the  solution.   Qualitative  as 
well  as  some  quantitative  aspects  of  the  solution  can  be  had 
through  a  study  of  the  locations  and  types  of  solution  curves 
existing  near  singular  points.   It  is  usually  desirable  to  have 
appropriate  equations  relating  the  variables  of  the  system  al- 
though  in  some  cases  relations  available  only  in  graphical  form 
may  be  used  (1,2). 

A  graphical  representation  of  solution  curves  on  a  plane 
surface  with  two  dimensions  is  conveniently  used  in  a  study  of 

singularities;  thus  such  a  study  is  limited  to  the  case  of  two 

dy   Q(x,y) 

variables.   If  a  differential  equation  of  the  form  —  =  

dx   P(x,y) 

is  investigated,  where  P ( x,  y)  and  Q(x,  y)  may  be  nonlinear 

functions  of  x  and  y,  the  equation  is  equivalent  to  the  two 

equations : 

dx  dy 

—  =  P(x,  y)         —  =  Q(x,  y) 

dt  dt 

and  is  obtained  by  eliminating  the  independent  variable  t. 

dy   Q(x,  y) 

Elimination  of  t  makes  —  =  an  autonomous  equation  and 

dx   P(x,  y) 

limits  it  to  situations  where  any  forcing  function  is  either 

entirely  absent  or  extremely  simple. 
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dy   Q(x,  y) 

Singularities  of  —  =  are  those  values  of  x  and  y 

dx   P(x,  y) 

for  which  both  P  and  Q,  become  simultaneously  zero.   At  the 

singular  point  x  =  xs  and  y  =  ys  and  P(xa,  ys)  =  0, 

dy 
Q(xg,  ys)  =  0,  and  —  becomes  indeterminate.   There  may  be  a 

dx 
number  of  singularities  if  P  and  Q,  are  nonlinear.   Since 

dx         dy 

—  =  0  and  —  =  0  at  a  singularity,  the  singularity  is  always  a 

dt         dt 

point  of  equilibrium.   As  was  noted  in  the  introduction,  the 
nature  of  solutions  near  a  singularity  may  be  explored  by  sub- 
stituting x  =  xs  +  u  and  y  =  ys  +  v,  where  u  and  v  are  small 
variations.   With  these  substitutions, 

dy   Q(x,  y) 


dx   P(x,  y) 
becomes : 

dy    dv   Q(xs,ys)  +  Cu  +  Dv  +  C2u2  +  D2v2  +  F2uv 


dx    au    P(xs,us)  +  Au  +  Bv  +  A2u2  +  B2v2  +  E2uv  + 


S>  ^S 

A,  B,  C,  D,  real  constants.   A  Taylor's  series  expansion 

may  be  necessary.   The  most  important  terms  in  determining  the 
solution  near  a  singularity  are  the  linear  terms  in  u  and  v. 
The  kind  of  singularity  depends  only  upon  the  linear  terms, 
provided  these  terms  are  present  if  nonlinear  terms  are  present. 
In  other  words,  if  C2  j   0  in  the  numerator  so  that  a  term  u2 
appears,  the  linear  term  Cu  must  be,  present  with  C  j   0  for  the 
singularity  to  be  simple.   The  same  condition  applies  in  both 
numerator  and  denominator  for  both  u  and  v.   Under  this  condition 
with  a  simple  singularity,  the  properties  of  the  solution  near  a 
singularity  depend  upon  the  equation: 


58 


dv    Cu  +  Dv 


(1) 


au    Au  +  Bv 
Only  linear  terms  appear.   If  there  are  higher  power  terms  in 
u  and  v  present,  a  study  cannot  be  made  from  this  equation  alone. 
Equation  1  is  equivalent  to  the  pair  of  simultaneous  first- 
order  equations: 

dv 


du 

—  =  Au  +  Bv 

dt 


=  Cu  +  Dv 


dt 


This  pair  of  equations  is  a  simple  case  of  the  more  general  set 
of  n  simultaneous  first-order  equations. 

A  slight  digression  here  for  a  discussion  of  a  technique 
which  simplifies  matters  later  is  in  order.   The  set  of  n  first- 
order  equations  may  be  written 

xl  =  allxl  +  a12x2  +  *  '  '  alnxn 


x2  =  a21xl  +  a22x2  +  '  •  •  a2nx 


n 


(2) 


*n  ~  anlxl  +  an2x2  +  •  ■  •  annxn 


where  x-,  . 
Moreover : 


.  x   are  the  n  dependent  variables 


a]    = 


x  = 


dx 

dt 


(x)      -    K} 


(3) 


'11      a12    *     ' 


.    a 


In 


'nl      an2 


•       •       •      8 


nn 


59 


Solutions  for  equations  of  the  form  (2)  are  well  known  and 
involve  exponential  functions  which  retain  their  form  upon  dif- 
ferentiation.  Thus 

fxj   =  {  cx]  exP  (^ 
A.  is  a  constant  determined  by  the  coefficients  in  the  differ- 
ential equations.   Differentiating  the  solution  and  substituting 
into  (3)  yields 

X(x)  =  |A]{x} 

which  may  be  written 

1 


[a]-  x   [i]J=  fo) 


J 

where  I  is  the  identity  matrix.   The  equation  may  be  satisfied, 

except  for  the  trivial  case,  only  if 

r 


lM 


x 


il 

L  J 


=     O 


where  the  A's  are  the  characteristic  roots  or  eigenvalues. 

Physical  systems  can  often  be  described  by  differential 
equations  of  the  form  of  (l)  and  usually  the  right  sides  of  the 
equations  involve  several  of  the  dependent  variables.   A  system 
in  which  there  is  no  coupling  would  have,  on  the  right  side  of 
each  equation  only  the  single  variable  which  appears  on  the 
left.   Such  an  equation  is  said  to  be  in  normal  form.   A  system 
with  coupling  may  be  converted  to  one  having  no  coupling  by  the 
mathematical  process  of  changing  the  variables  through  an  ap- 
propriate linear  transformation.   Thus 

xl  =  Pil*l  +  Pl2^2  +  •  •  .  Pln^n 
x2  =  P21S1  +  P22^2  +  •  '  •  P2n^n 


*n  =  Pnl^l  +  Pn2^2  + 


Pnn-^n 
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In  matrix  form,  |xj  =  j_BJ  /y|  ,  P-j_  ^  are  constant  quantities. 
Rewritten : 


r-1 


{t)-  =  DT    W>   *l*° 


Rewritten : 


Thus: 


[pj^i)    =    {yj    =  pj1   [AJ  (x)     =  QP]-1  pQ  ^    (yj 


(U) 


A  solution  for  (I4.)  has  the  form  jy|  =  \CYJr    exp(At)  and 
the  characteristic  roots  of  (I4.)  are  given  by:   [j3J  -  X«  rij  =  0. 


Thus 


<-l 


,-1 


[B]  -  ,B  [I]  =  [p]   [A]  [P]  -  ,B  |V|   [I]  [P] 

=  w1  [m  -  ^  pa]  m 

The  determinant  of  both  sides  must  vanish.   Since   P   f  0  and 
P"1   f  0,  it  must  be  that  [jT\   -  AB  |~I^[  =  0.   Therefore  the 
transformation  [_Pj  changing  [~-A~j  to  [~B"j  does  not  change  the  set 
of  roots  for  the  set  of  equations. 

"If  the  transformation  matrix  PPJ  is  chosen  properly,  [_B_{ 
can  be  made  a  diagonal  matrix  with  all  the  elements  not  on  the 
main  diagonal  zero.   The  set  of  equations  described  by  QbJ  are 
then  said  to  be  in  normal  form.   The  eigenvalues  of  I  Bj  are  the 
same  as  those  of  FaJ  .   The  elements  on  the  main  diagonal  of  [j3j 
must  be  the  eigenvalues.   Thus  the  set  of  equations  in  normal 
form  is: 

71  =  *i7i 

y2  =  \2J2 


yn  =  ^n7n 
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Lb]  = 


X- 


0 


0 
X, 


0 
0 


0 


A 


n 


If  two  roots  happen  to  form  a  complex  pair  as   X-^  =  6  +  jco, 
X2  =  6  -  jco,  these  two  elements  of  |_B  J  would  be  complex  quanta 
ties.   Often  it  is  desirable  to  have  only  real  quantities  ap- 
pear.  An  equivalent  form  for  the  equation  when  X-j_  and  X2  are 
complex  quantities  is 

* 

71  =  12 

y2  =  -(52  +  a)2)y1  +  2  5y2 


[b]  = 


Jn  =  Xn? 


n 


0  1 

-(62  +  go2)    25 


0 


0 
0 


X 


n 


The  alternate  form  is  usually  preferable  when  complex  roots 
occur.   The  original  set  of  equations  jxj  =  lAj/xj  ,  has  solu- 
tions {xj  =  jCxl  exp(Xt)  .   In  normal  form  f  yj  =  (Cyj  exp(Xt). 
The  only  feature  of  the  solutions  which  can  be  determined  from 
the  differential  equations  alone  are  the  eigenvalues  which  are 
the  same  in  each  case.   The  coefficients  (Cxf  and  jCyl  can  be 
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found  only  from  initial  conditions.   Solution  in  terms  of  x  and 
y  can  then  be  said  to  be  equivalent,  at  least  so  far  as  quali- 
tative properties  are  concerned. 

Now  the  second-order  system  is  described  by  just  two  equa- 
tions.  The  matrix  [_PJ  needed  to  reduce  this  kind  of  system  to 
normal  form  is  quickly  found.   Consider  the  use  of  the  follow- 
ing symbols. 


W-O  m: 


b 

d 


t» 


*    (3 
6 


(y? 


i 


■*1 

y2j 


[B]  = 


A- 


0 


0 


\< 


m  -  *  dh 


=  0 


Originally,  only  |x)  and  [~Aj  are  known,  and  [_Bj  ,  [_PJ  ,  and  fyl 
must  be  found.   The  characteristic  equation  is 

a  -  A     b 
c     d  -  A 

or   A^    (a  +  d)A  -  (be    ad)  =  0.   The  two  characteristic  roots 
are : 

1  (  r       2  -i  1I2 
(A1,  A2)  =  -   (a  +  d)  ±1  (a  +  d)   +  ij.(bc  -  ad) 

2  ( 

where,  in  all  that  follows,  A-,  is  the  root  found  with  the  posi- 
tive sign.   If  the  roots  are  real,  A-,  is  always  the  more  positive 
root.   These  values  of  X-,  and  A2  serve  to  determine  matrix  [j3J  . 
Matrix  [V]  must  satisfy  [bJ   =  []PJ  [Y]  [V]  .   Therefore 

|_P]]  JJETJ  =  [~A^j  £p^  .   When  these  matrix  products  are  found,  the 
result  is: 
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*x- 


i\- 


p\2 

6X2 


-<a  +  Vb         pa  +  6b 
-<c  +  Vd    pc  +  6d 


Since  corresponding  elements  of  equal  matrices  must  be  iden- 
tical, the  following  simultaneous  equations  exist: 

/ 

-<\]_  =  -<a  +  /b  ;  Y^i   ~   "<c  +  ^d'   The  ratio  —   defined  as  m-^  is: 

//   X-i  -  a     c 


mi  =  —  = 


+  b 

5 

P 
X2  ~    a 


A.-L  -  d 


In  a  similar  manner,  the  ratio  —  is  defined  as  m2  and  is: 


mo  =  — 


A2  -  d 


p     b 

Both  forms  are  necessary  on  the  right-hand  side  of  the  above 

equations  since  sometimes  one  form  is  indeterminate.   The  ratios 

y  6 

—  and  —  are  fixed  by  these  equations  and  thereby  fix  the  ele- 

=<      p 

ments  of  matrix  P  within  constant  factors.   The  coordinate 

transformation  is  fxj  =  j_PJ  jyj  ,  or 
*1  =  ^1\   +    Py2 

*2  =  ?7l   +  &J2   =  ^l^Jl   +  ^2^2 
The  types  of  singularities  of  a  second-order  system  can  now 

be  investigated  and  classified.   -The  simplest  cases  are  those  in 

which  the  two  characteristic  roots  are  real.   Equations  for  the 

system  in  normal  form  are: 


0 


0 


\. 


and 


6k 


y2  •   ^2^2    d^2 

yi     xiyi     d^i 

Possibilities  may  be  listed  as  follows  where  X-,    designates 

the  more  positive  root. 

1.  Both  Roots  Real  and  Positive . 

A2 
0   <  X2  <   X1    ;  0  <  —  <  +  1 

The  equation  for  a  curve  representing  a  solution  for  the  dif- 
ferential equation  on  the  y-,yp  plane  can  be  found  directly  by 
integration: 

C  is  an  arbitrary  constant  dependent  upon  initial  conditions. 
These  curves  are  generally  parabolic  in  shape,  with  the  exact 

x2 

shape  determined  by  —  and  constant  C.   The  slope  of  the  curves 

dy2      ^2    ((xP  \n)-l)  .  . 

may  be  found  from  =  C  —  y-,    c      ±  and  near  the  origin, 

dy]_      ^! 

dy2  A2 

?•  oo   as  j-\ *  0,  since  —  <  +  1.   Thus  all  solution 

dyi  Xx 

curves  have  a  definite  direction  near  the  origin,  being  paral- 
lel to  the  y2  axis.   Shown  below  is  the  case  where  values  of 

\2 
A.  2  and  A.2  are  not  far  different,  i.e.,  0  << —  <1.   The  curves 

Xl 
represent  the  locus  of  points  determined  by  corresponding 

values  of  y-^  and  y2.   As  independent  variable  t  increases,  the 
point  relating  instantaneous  values  of  y,  and  y2  moves  along 
the  curve  in  the  direction  of  the  arrowheads.   Initial  condi- 
tions determine  the  value  of  constant  C,  and  thus  the  quadrant 
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-+■ 


•  y- 


within  which  a  particular ' solution  lies.   Since  the  roots  are 
positive,  both  y-,  and  yp  increase  without  bound  as  t  increases 
and  this  type  of  singularity  is  said  to  be  unstable.   On  the 
other  hand,  if  y-,  and  yp  were  both  to  vanish  as  t  increased,  the 
singularity  would  be  called  stable.   One  of  the  important  fea- 
tures of  a  singular  point  is  the  question  of  its  stability. 
Node  is  the  name  given  to  this  type  of  singularity,  referring  to 
the  fact  that  the-  solution  curves  have  a  definite  direction  near 
the  singularity.   Either  the  y-,  or  y?  axis  can  be  a  solution  if 
the  initial  conditions  are  such  that  one  of  the  variables  y-j_  or 
yp  is  exactly  zero.   The  axes  then  represent  special  cases  of 

solution  curves  to  correspond  to  special  initial  conditions. 

X2 
If  the  values  of  \-^_   and  Ap  are  such  that  0  <  —  <r<l,  "the 

xl 
solution  curves  take  the  shape  of  the  figure  below.   These  curves 

show  a  similarity  to  the  curves  above,  but  as  they  change  direc- 
tion the  breaks  are  much  sharper. 


A  ^2 


t 


■*•  y 


i 


66 


Consider  the  following  solutions  of  the  form: 

y1  =  Cj_  exp(X1t)  ;       y2  =  C2  expU2t) 

dy2    A2C2  exp(A.2t) 

=  0  <  X2  <  <  A.]_ 

dy-j_   ^ici  exp(X1t) 

If  t  is  large  and  positive: 

\it  .  w  Xpt    ,   ^2 
e  -1-  >  >  e  d      and  &  0 

dyx 
If  t  is  large  and  negative: 

A. it  ^    j-    \ot         ,   y2  ^ 
e  x  <  <  e  *■      and  —  P° 

Thus  if  X-j_  and  A2  differ  sufficiently  in  magnitude,  transition 
between  these  two  occurs  suddenly  and  solution  curves  are  essen- 
tially two  straight  lines  with  a  sharp  break  joining  them. 
2.  Both  Roots  Real  and  Negative . 

x2 

\n  «=  \i  '<  0   ;      —  >   +  1 


■1 


*1 


Solution  curves  for  the  normal  form  are  again  parabolic, 

dy2 
but  near  the  origin  =  0  and  the  curves  are  parallel  to  the 

dy; 
y-.  axis.   The  singularity  is  again  a  node,  but  since  negative 

eigenvalues  lead  to  ultimate  disappearance  of  y-,  and  y2  with 

increasing  t,  the  node  is  stable.   For  increasing  t  the  point 

representing  corresponding  values  of  y-,  and  y?  moves  closer  to 

the  singularity  at  the  origin  but  reaches  it  only  at  the  infinite 

value  of  t.   The  y-,  and  y?  axes  are  again  special  solution  curves 

3.  Both  Roots  Real  and  Opposite  in  Sign. 

x2 

X?   <  0  <  X-,  —  <  0 
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In  this  case 


dy2 


X2    12  U?Ai  I 

—   ( — )  and  y2^l        =  c-   Solution 


curves  in  the  y^y2  plane  are  hyperbolic  in  shape  and  generally 
pass  by  the  singularity  at  the  origin.   Now  X-^_   i-3   positive  and 
y,  ultimately  increases  without  bound  and  the  solution  is  un- 
stable, even  though  y^  ultimately  vanishes.   This  singularity 
is  the  saddle. 

Lj..  Roots  Pure  Imaginaries . 

A-l  =  +jco 

\2   ~   -J'w 


5  =  0 


d^2      9  Yl 
The  equation  for  a  solution  curve  can  be  found  from  =  -co^  — 

dyi       y2 

?      o  o 

as  co  y-i   +  Jn      =  C;  the  equation  is  that  of  an  ellipse  about  the 

singularity  at  the  origin.   The  figure  below  shows  a  general 
case.   This  type  of  singularity  is  known  as  a  vortex.   The  solu- 
tion is  a  periodic  oscillation  in  time  with  no  change  in  amplitude 

4  ^2 


*1 


There  is  neither  growth  nor  decay;  the  solution  has  a  "neutral" 
stability  and  the  amplitude  and  size  of  the  ellipse  are  deter- 
mined by  the  initial  conditions. 

5-  Roots  Complex  Conjugate  s. 

A-[_  =  6  +  jo)  ;       A2  =  6  -  jco 

In  this  case  the  normal  form  for  the  equations  can  be 
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written 


0       1 
-(S2  +  GO2) 


1 
26 


yi 


dy2    -(62  +  0)2)y1  +  2  5y2 


dyx 


y 


and  solution  curves  depend  upon  this  equation.   Integration  is 
possible  following  an  appropriate  change  of  variable.   The  qual- 
itative nature  of  its  solution  curves,  however,  may  be  obtained 
more  easily  from  an  observation  of  lines  of  equal  slope.   The 

following  seem  evident  from  the  equation: 

2' 


dy2 
dyi 
dy2 
dyi 
dy2 
dyx 


=  0    along   y2  = 


(5^  +  oo^)y1 


25 


-  oC      along    y2  =  0 


=  25   alon 
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These  isoclines  carrying  directed  line  segments  of  appropriate 
slope  are  shown  below  and  typical  solution  curves  sketched. 


4    y? 


yi 


3>o 


§<:o 
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The  curves  form  spirals  about  the  singularity  at  the  origin 
which  is  designated  as  a  focus.   If  6   0,  the  solution  ulti- 
mately grows  without  bound  and  is  unstable.   If  6   0,  the 
solution  ultimately  vanishes  and  is  stable. 

For  simple  singularities,  then,  of  two  first-order  linear 
equations,  there  are  only  four  possibilities:   node,  saddle, 
vortex,  and  focus.   The  node  and  focus  may  be  either  stable  or 
unstable,  the  saddle  is  always  unstable,  and  the  vortex  is 
"neutrally"  stable.   Only  the  types  of  solution  curves  associ- 
ated with  these  four  singularities  can  exist. 
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Three  approaches  to  investigations  of  nonlinear  stability 
are  described  in  this  report.   The  three  approaches  are:  The 
Describing  Function  Method,  Analysis  by  Means  of  Singular 
Points ,  and  The  Variable  Gradient  Method  of  Generating  "V" 
Functions  used  in  conjunction  with  certain  theorems  attribut- 
able to  Lyapunov. 

In  Appendices  A  and  B  some  derivations  and  descriptions 
are  contained  which  provide  background  material  applicable  to 
the  describing  function  and  singularity  point  analysis  ap- 
proaches, respectively.   The  main  body  of  the  report  was  then 
utilized  to  exhibit  some  results  of  the  various  approaches 
applied  to  nonlinear  stability  problems. 


